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ABSTRACT 

Numerical  solution  to  a  theoretical  model  of  vapor  cavitation  in  a 
dynamically  loaded  journal  bearing  is  developed,  utilizing  a  multigrid 
iterative  technique.  The  method  is  compared  with  a  noniterative  approach  in 
terms  of  computational  time  and  accuracy.  The  computational  model  is  based  on 
the  Elrod  algorithm,  a  control  volume  approach  to  the  Reynolds  equation  which 
mimics  the  JakoOsson-Floberg  and  Olsson  cavitation  theory.  Besides  accounting 
for  a  moving  cavitation  boundary  and  conservation  of  mass  at  the  boundary,  it 
also  conserves  mass  within  the  cavltated  region  via  a  smeared  mass  or  striated 
flow  extending  to  both  surfaces  in  the  film  gap.  The  mixed  nature  of  the 
equations  (parabolic  in  the  full  film  zone  and  hyperbolic  in  the  cavitated 
zone)  coupled  with  the  dynamic  aspects  of  the  problem  create  interesting 
difficulties  for  the  present  solution  approach.  Emphasis  is  placed  on  the 
methods  found  to  eliminate  solution  instabilities.  Excellent  results  are 
obtained  for  both  accuracy  and  reduction  of  computational  time. 

NOMENCLATURE 

AR  aspect  ratio  of  grid  size,  Ax/Az 

e<j  dynamic  eccentricity,  m 


static  eccentricity,  m 
forcing  function  on  grid  k 
residual  function  on  grid  k 
switching  function 
dimensionless  film  thickness,  h/AR 
f i lm  thickness,  m 

interpolation  function  from  grid  k  to  grid  k-1 
grid  indicator,  k  <  M 

differencing  scheme  acting  on  a  variable  6 

1  ength-to-dl ameter  ratio 

represents  the  finest  grid  (highest  number) 

number  of  grid  points  axially 

1  i neal  mass  flux ,  kg/m-s 

number  of  grid  points  circumferentially 

fluid  pressure,  N/m2 

ambient  pressure,  N/m2 

cavitation  pressure,  N/m2 

radius  of  journal ,  m 

time,  s 

sum  of  the  surface  velocities  in  x-direction,  m/s 
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AR  radial  clearance,  m 

At  time  increment,  s 

Ax  incremental  spacing  along  circumference,  m 

Az  axial  Incremental  spacing,  m 

c  eccentricity  ratio 

0  fractional  film  content  In  cavitated  zone 

density  ratio,  p/pc,  in  full  film  zone 

p  dynamic  viscosity,  N-s/m^ 

p  fluid  density,  kg/m^ 

pc  fluid  density  within  cavitated  zone,  kg/m^ 

angular  coordinate  along  circumference,  degree 

wd  orbital  angular  velocity  of  journal  center  about  a  fixed  point 
relative  to  the  housing  center,  rad/s 

015  angular  velocity  of  journal  about  its  own  center,  rad/s 

INTRODUCTION 

The  presence  of  vapor  cavitation  in  dynamically  loaded  journal  bearings 
has  become  a  topic  of  increasing  Importance.  The  use  of  increased  loads  and 
more  complicated  loading  cycles  has  resulted  in  an  increase  in  the  occurrence 
of  cavitation  erosion  problems.  Examples  of  journal  bearing  applications 
include  main  and  crankshaft  bearings  in  diesel  engines  and  a  variety  of 
bearings  In  the  aircraft  industry,  Dowson  and  Taylor  (]_)•  Dynamic  loading  can 
also  lead  to  instabilities  in  the  motion,  such  as  whirling  or  whipping  motion, 
which  may  damage  the  bearing.  In  order  to  avoid  bearing  damage,  it  is  useful 
to  predict  the  conditions  under  which  the  bearing  will  remain  stable.  The 
determination  of  these  stability  maps  requires  a  knowledge  of  the  hydrodynamic 
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Previous  static  loading  models,  sued  as  those  using  Swift-Stieber  or 
Giimbel  boundary  conditions,  assume  a  stationary  cavitation  bubble  and  are 
inadequate  for  high  speed,  dynamic  applications.  Under  dynamic  loading, 
changes  in  the  local  film  thickness  cause  the  bubble  to  grow,  move  downstream 
from  the  minimum  film  position,  and  collapse,  Brewe  (2). 

A  film  model  which  effectively  deals  with  dynamic  loading  has  been 
formulated  by  Jakobsson-Floberg  (3)  and  Olsson  (4).  Besides  accounting  for  a 
moving  boundary,  it  also  accommodates  the  flow  within  the  cavitated  region, 
manifested  via  a  smeared  mass  or  liquid  striations.  These  striations  have  been 
observed  in  past  experimental  work.  Because  the  theory  assumes  a  zero  pressure 
gradient  within  the  cavitated  zone,  this  mass  flow  is  a  Couette  flow.  The  JFO 
theory  accounts  for  both  film  rupture  and  film  reformation,  another  advantage 
over  previous  methods.  Unfortunately,  the  complexity  of  the  JFO  theory  makes 
it  difficult  to  apply  (2).  Elrod  and  Adams  (5,6)  have  developed  an  algorithm 
which  automatically  conforms  to  the  JFO  theory  while  being  much  simpler  to 
code.  It  utilizes  a  switching  function  which  eliminates  the  pressure  gradient 
terms  from  the  Reynold's  lubrication  equation  at  cavitated  points. 

A  solution  to  the  Elrod  algorithm  for  a  dynamically  loaded  problem  has 
been  formulated  by  Brewe  (2).  This  direct  solution  to  the  finite  differenced 
equations,  i.e.,  no  iterations  required,  utilizes  an  alternating  direction 
implicit  (ADI)  scheme  for  the  time  march.  When  compared  to  a  nonconservative 
film  model  (pseudo-Giimbel  boundary  conditions),  as  much  as  a  20-percent 
difference  in  load  capacity  is  observed.  Brewe's  results  agree  excellently 
with  the  experimental  work  of  Jakobsson  and  Floberg  (3)  for  stationary 
cavitation.  The  present  experimental  data  on  nonstationary  cavitation  is 
limited,  but  Brewe's  results  compare  reasonably  well  with  the  experimental 
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work,  of  Jacobson  and  Hamrock,  (7,8).  Unfortunately,  Brewe's  direct  method 
requires  two  to  three  times  the  computational  work  needed  by  the  nonconservative 
solution.  The  practical  use  of  the  Elrod  algorithm  in  solving  dynamic  loading 
problems  in  industry  requires  the  development  of  a  more  efficient  computer 
solution. 

Iterative  techniques  are  not  considered  among  the  fastest  methods  of 
computer  solution.  However,  newly  developed  techniques  using  multiple  grid 
sizes  have  shown  that  It  is  possible  to  greatly  reduce  computational  time.  It 
is  the  purpose  of  this  work  to  Implement  the  multigrid  technique  developed  by 
Achi  Brandt,  (9,J_0) ,  in  the  iterative  solution  of  the  Elrod  algorithm  under 
dynamic  loading  conditions.  As  a  point  of  reference,  the  multigrid  method  has 
also  recently  been  implemented  In  an  EHD  lubrication  problem  by  Lubrecht, 
ten  Napel  ,  and  Bosma,  ( 1J_) . 

The  nature  of  the  problem,  i.e.,  the  presence  of  discontinuous 
coefficients  In  the  cavitated  region,  poses  interesting  difficulties  in  the 
application  of  the  multigrid  method.  Therefore,  one  major  objective  is  to  find 
methods  of  making  the  multigrid  technique  as  effective  over  the  cavitated  area 
as  it  is  over  areas  of  full  film. 

DESCRIPTION  OF  THE  MODEL  PROBLEM 

The  bearing  motion  consists  of  a  journal  undergoing  a  constant  rpm  as  well 
as  a  noncentered  circular  whirl  inside  of  a  360°  cylindrical  bearing.  The 
journal  center  moves  through  a  prescribed  dynamic  cycle,  (Fig.  1)  from  a 
minimum  through  a  maximum  eccentricity  (see  Table  I  for  operating  conditions). 

The  theoretical  model  assumes  conditions  of  heavy  loading,  i.e.,  load 
carrying  capacity  >>  surface  tension  forces  in  the  liquid.  An  oil  lubricant 
is  used  for  which  the  vapor  pressure  is  very  nearly  zero.  It  is  assumed  that 
the  lubricant  has  been  degassed  and  that  only  vapor  cavitation  is  ptesent.  In 
actual  experimental  work  with  submerged  bearings,  this  is  a  necessary 


condition,  since  gaseous  cavitation  forms  at  near  ambient  pressure  and  might 
prevent  the  occurrence  of  subambient  pressures.  The  flow  balance  into  and  out 
of  the  bearing  cannot  be  maintained  in  the  absence  of  subambient  pressures. 
Within  the  cavitated  region,  a  zero  pressure  gradient  is  assumed.  In  order  to 
determine  the  load  carrying  capacity  of  the  bearing  under  the  prescribed 
motion,  the  flow  equations  are  solved  for  the  pressure  variable.  It  is  more 
convenient,  however,  to  introduce  another  dependent  variable,  9,  which  has  a 
dual  interpretation  regarding  the  full  film  region  and  the  cavitated  region. 

In  the  full  film  (9  2  1-0),  9  Is  the  ratio  p/pc.  which  represents  the  ratio 
of  film  mass  content  to  the  mass  content  that  would  exist  at  the  cavitation 
pressure  pc.  In  the  cavitated  region  (9  <  1.0),  9  represents  the  mass 
content  that  exists  in  the  form  of  liquid  striations.  Pressure  and  density 
are  related  through  the  bulk  modulus,  3,  such  that  (p3p/8p)  =3-  A  switching 
function  (g(6) )  automatically  eliminates  the  pressure  term  at  cavitated  points 
of  the  flow.  That  is, 

Uncavitated  point  (9  >  1.0):  g  =  1 

Cavitated  point  (9  <  1.0):  g  =  0 

THE  GOVERNING  EQUATION  AND  DIFFERENCING  SCHEME 

The  Reynold's  lubrication  equation  written  In  terms  of  the  fractional  film 
content,  the  switch  function,  and  the  bulk  modulus  becomes 

*  (!)  •  ♦<»>  -  ?  *  (?£)«>* 

Note  that  the  right-hand  side,  the  pressure  induced  Flow,  completely 
disappears  in  the  cavitated  region  where  the  switch  function  becomes  zero.  The 
finite  difference  equations  are  obtained  using  a  control  volume  approach  as  was 
used  by  8rewe  (2).  A  control  volume  (Fig.  2)  is  constructed  about  each  nodal 
point  and  the  net  change  in  mass  flow  into  the  cell  (mx+^x/2  -  mx-Ax/2^  is 
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equated  to  the  total  Increase  of  mass  (pc9h)  In  the  cell  over  a  time  At.  The 
control  volume  equivalent  to  the  mass  conservation  equation  is 


Am  Am 
_ x  _ z 

Ax  +  A2 


In  the  full  film  region  (all  g  =  1),  both  the  convective  and  pressure 
terms  are  central  differenced,  appropriate  for  the  parabolic  system.  In  the 
fully  cavitated  region  (all  g  =  0),  the  pressure  terms  are  eliminated  and  the 
convective  terms  are  upwind  differenced,  9  accounting  for  the  mass  transport. 
The  combination  of  switching  terms  at  the  cavitation  boundary  automatically 
sets  well  posed  boundary  conditions  between  the  two  systems.  The  resulting  Am 
terms  from  the  control  volume  analysis  are 

(“x)conv  ■  "c©  Ih-l('  -  -  ho(>  -  9o)9o  *  5  klh-l(2  -  90) 

*  90ho(9-l  -  2  »  9,)  -  91l',90]j 

Wpress  ■  (r|)(k)  Sh-l/2  9-l(e-l  -  ')  -  (h-l/2  *  h?/2>o(eO  -  ') 

*  [h?/2  9,(e,  -  |)]( 

An  Euler  implicit  time  differencing  scheme  is  used  for  stability  purposes. 


press 


h-,/2 


giving: 


eh  -  e*h* 

At 

where  0*h*  signifies  time  t  -  1. 


Pc  Ax  pc  Az 


It  should  be  noted  that  all  terms  on  the  right-hand  side  of  this  equation 
are  evaluated  at  time  t  and  are  therefore  unknown.  Along  the  axial 
doundaries,  i.e.,  along  the  edge  of  the  bearing,  the  boundary  condition  is 
that  of  atmospheric  pressure.  The  circumferential  direction  has  wrap  around 
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boundary  conditions.  The  problem  and  boundary  conditions  are  described  In  more 
detai 1  in  (2) . 

THE  MULTIGRID  METHOD 

Analysis  of  the  Single  Grid  and  the  Residual  Function 

The  iterative  solution  of  a  set  of  equations  on  a  single  grid  generally 
has  rapid  convergence  over  the  first  few  sweeps,  but  very  slow  convergence  over 
most  of  the  process.  By  examining  the  solution  process,  the  reasons  for  thts 
become  clear.  Assume  the  following  continuous  differential  equation, 

L {©( x ) }  =  F(x);  with  suitable  boundary  conditions 
for  which  the  discretized  set  of  equations  on  one  grid  take  the  form 

Lk{6k}  =  Fk  (i) 

In  the  present  notation,  k  represents  the  particular  grid  size  used, 
represents  the  differencing  scheme  acting  on  9,  and  9  is  the  exact  solution 
of  the  differenced  equations. 

Let  9  be  the  present  approximation  to  the  exact  solution  9.  By 
substituting  9  Into  the  differenced  equations,  the  following  is  obtained. 

fk  =  Fk  _  Lk{9k}  (2) 

where  f^  is  referred  to  as  the  "residual  function."  The  residual  function 
is  a  means  of  analyzing  the  error  left  in  the  present  approximation. 

A  Fourier  analysis  of  fk  breaks  the  error  into  its  high-  and  low- 
frequency  components.  High  frequency  is  defined  as  wavelength  less  than  or 
equal  to  four  times  the  grid  spacing.  After  a  few  relaxation  sweeps,  e.g., 
Gauss-Sei del ,  these  high-frequency  components  are  smoothed  out,  due  to  the  fact 
that  they  are  locally  corrected.  Once  the  error  is  all  low  frequency,  the 
smoothing  rate  drops  drastically.  The  grid  spacing  is  too  fine  to  efficiently 
smooth  these  low-frequency  terms  Brandt  (9). 
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The  Roles  of  High-  and  Low-Frequency  Error  Smoothing 

The  basic  thrust  of  multlgrid  Is  to  utilize  coarser  grids  to  handle  these 
low-frequency  error  terms.  As  long  as  the  residual  function,  i.e.,  the  error, 
is  well  represented  on  a  particular  grid,  that  grid  can  quickly  smooth  its  own 
"high-frequency"  terms  and  send  an  appropriate  correction  back  to  the  finer 
grid.  By  utilizing  coarser  and  coarser  grios,  all  of  the  low  frequency  terms 
are  treated  similarly.  Besides  their  ability  to  efficiently  deal  with  low- 
frequency  terms,  coarse  grids  also  have  fewer  nodal  points  to  sweep  through, 
making  the  coarse  grid  sweeps  very  cheap. 

Whereas  moving  to  coarser  grids  smooths  the  low-frequency  error,  fine  grid 
updates  during  the  multigrid  process  have  value  in  improving  the  accuracy  of 
the  present  solution.  The  role  of  relaxation  on  the  finer  grid  is  to  resolve 
its  own  high-frequency  components  as  well  as  smooth  the  high-frequency  error 
which  is  produced  by  interpolation  from  the  coarser  grid. 

Coarse  Grid  Representation 

The  fine  grid  problem  itself,  i.e.,  =  F^,  is  not  what  is  really 

being  represented  on  the  coarse  grid  k-1.  The  actual  purpose  of  the  coarse 

k-1 

grid  is  to  solve  for  a  correction  value,  0C 

e£_1  =  ek_!  -  9k~' 

as  a  function  of  the  amount  of  residual  .rror  that  is  left  in  the  present 

approximation  on  the  fine  grid,  i.e.,  fk. 

For  a  nonlinear  problem,  the  full  approximation  storage  (FAS)  mode  must  be 

~k-l 

used.  This  method  stores  the  entire  value  of  9  on  the  coarse  grid  k-1 

k- 1 

instead  of  just  the  correction  value  0C  .  If  used  on  a  linear  problem,  the 
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equations  reduce  to  those  for  the  linear  mode  Brandt  (9).  The  existence  of 
nonl Inearl  ties  In  the  convective  mass  flow  contribution  necessitates  the  use 


of  the  FAS  mode.  Using  FAS,  the  coarse  grid  problem  becomes 

rk-1  t k-l/sk 

initial  "  lk  19 


where  I 


Also, 


Solve 


k-l 


is  the  interpolation  operator  from  the  fine  to  the  coarse  grid. 


=  I 


'(fk  -  ikek) 


+  Lk-1{ek_1} 


Lk-1{ek_1}  =  Fk_1 


(3) 


e 


k-l 


ek_1 


e 


k-l 

Initial 


Using  the  interpolation  operator  lj^_p  i.e.,  from  the  coarse  to  the  fine 
grid,  we  obtain 

■k  ,(ek-’) 


*  §otd  *  ’k-r 


(4) 


k-l  k 

The  fact  that  the  coarse  grid  is  solving  for  9c  as  a  function  of  f 

ir 

has  two  Important  consequences.  The  first  is  that  f  must  be  well  represented 

k-l 

on  the  coarse  grid  in  order  that  9  is  an  accurate  correction  value. 

c 

Jr 

Therefore,  care  must  be  taken  that  f  Is  well  smoothed  on  the  fine  grid  before 
transferring  it  to  the  coarse  grid. 

The  second  consequence  Is  that  It  is  not  necessary  that  the  coarse  grid 
k-l  ~k-l 

differencing  scheme,  L  {9  },  exactly  match  the  fine  grid  differencing 

K  ~k 

scheme,  L  {9  } .  This  can  be  seen  from  the  following  consideration,  i.e.,  the 


:oncept  of  solving  for  a  correction  value  on  the  coarse  grid  arises  because 

>] 


lLk-'.0k-'i  -  it"'#-1 


should 

represent 


Lk{9k)  -  Lk{9k} 


It  can  be  shown  that  the  multigrid  equations  follow  from  this  basic 
concept.  Now  considering  the  RHS,  ©k  is  unknown,  although  we  can  use  the 
fact  that  Lk{©k}  E  fk.  Therefore, 


RHS  =  [pk  -  lk{©k 


{©k}  = 


f  =  Residual  function. 


~k-l  k-1 

From  the  definition  of  ©c  ,  and  the  fact  that  L  is  a  linear  operator, 
the  LHS  can  be  wr i tten 

Lk_1{ek-1}  =  Lk_1 {ek_^ }  -  Lk_1{eK_1} 

The  full  equation  then  becomes  the  standard  linear  multigrid  equation 


k-1  r/ik- 1  ■ 


L  }  =  I* 


k-1 /A 


which  can  then  be  adapted  to  the  nonlinear  form  expressed  in  Eq.  (3). 

Therefore,  the  terms  which  must  closely  represent  each  other  are  the 

pracketed  terms,  not  the  individual  components  within.  Thus  there  is  some 

v_1  ~k_i 

flexibility  in  creating  l  {©  }.  Especially  if  the  problem  contains 

rapidly  changing  spatial  coefficients,  the  coarse  grid  differencing  scheme  will 
nave  to  be  a  modified  version  of  the  fine  grid  scheme. 

Mul ti grid  Cycle 

Various  multigrid  cycles  can  be  used.  When  developing  a  multigrid  code, 
it  is  best  to  use  a  prescribed  cycle  so  that  the  results  obtained  by  testing 
different  relaxation  and  interpolation  schemes  can  be  easily  compared.  One 
example  of  this  type  is  the  V-cyde.  One  V-cycle  consists  of  the  following: 
a  predefined  number  c  sweeps  on  each  grid  in  descending  order  of  fineness 
until  the  coarsest  grid  is  reached;  iterating  the  coarsest  grid  problem  to 
convergence;  and  a  predefined  number  of  sweeps  on  each  grid  in  ascending  order. 
The  goal  in  multigrid  is  to  obtain  an  order  of  magnitude  error  reduction  per 
cycle.  Once  tne  pest  relaxation  and  interpolation  schemes  for  the  problem  have 


been  determined,  the  adaptive  multigrid  method  can  be  used.  This  Is  generally 
more  efficient  than  the  prescribed  cycle  method.  The  adaptive  algorithm  stays 
on  a  particular  grid  until  convergence  on  that  grid  has  slowed  to  a  defined 
"slow"  rate,  at  which  time  it  automatically  moves  to  the  next  coarser  grid. 
Whenever  it  converges  to  the  set  tolerance  at  a  particular  grid,  it  moves  back, 
up  to  the  next  finer  grid. 


Determination  of  Smoothness  Between  Grids 


The  process  of  determining  smoothness  between  grids  is  described  using  FAS 
and  the  adaptive  multigrid  algorithm.  Let  M  denote  the  finest  grid. 


lV)  =  fm 


After  each  relaxation  on  grid  M,  the  program  decides  whether  to  stay  on 
the  present  grid  or  move  to  a  coarser  one.  For  simplicity  sake,  this  is  often 
done  by  measuring  the  rate  of  convergence  and  determining  a  cutoff  rate,  i.e., 
if  convergence  is  "slow,"  the  program  will  move  to  a  coarser  grid  (see 
Fig.  3(a)  for  a  schematic  of  the  general  problem).  This  method  Infers  that 
slow  convergence  signifies  a  smoothed  residual  function.  It  assumes  that  the 
presence  of  high-frequency  terms  will  show  up  in  a  rapidly  decreasing  global 
error  term,  where 


global  error  = 


-  6,  , 


ol  d 


i  ,J 


new. 


This  is  not  a  bad  assumption  if  there  are  no  local  areas  containing 
-apidly  changing  spatial  coefficients.  Problems  may  occur  if  local  regions  of 
nigh-frequency  residuals  exist  within  a  globally  smooth  domain.  The  global 
error  term  may  not  be  affected  by  the  high-frequency  errors  and  will 
interpolate  the  solution  to  the  coarser  grid,  where  the  residual  error  of  the 
local  region  will  not  be  well  represented.  Some  other  method  of  determining 
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the  smoothness  of  this  local  area  is  needed  in  such  a  case.  Additional  sweeps 
over  the  fine  grid,  or  parts  of  the  fine  grid,  are  then  needed  to  smooth  the 
error  locally. 

Assuming  that  the  residuals  have  been  sufficiently  smoothed  on  the  fine 
grid  k,  the  problem  moves  to  grid  k-1  and  relaxes  Eq.  (3).  The  program 
either  returns  to  the  fine  grid  if  the  solution  is  converged,  or  goes  to  a 
coarser  grid  if  the  convergence  is  slow  and  the  residuals  are  smooth.  When 
the  coarsest  grid  is  reached,  a  converged  solution  is  obtained  by  continued 
relaxation  or  by  direct  solution.  When  a  converged  solution  is  obtained  on  a 
grid  k-1,  the  solution  is  updated  on  grid  k  using  Eq.  (4). 

APPLICATION  OF  MULTIGRID 

The  use  of  Implicit  time  differencing  necessitates  the  solution  of  an 
NlxMJ  set  of  equations  at  each  time  step.  In  this  study,  no  attempt  is  made 
to  use  multigrid  across  physical  time.  Multigrid  is  used  to  facilitate  the 
iterative  solution  of  the  NIxMJ  set  of  equations  within  each  time  step. 

The  fine  grid  (M)  equations  take  the  form  LM{9M}  =  FM,  where  LM{6M} 
represents  the  differencing  scheme  of  9  described  earlier.  The  forcing 
function,  FM,  represents  terms  from  the  previous  time  step  which  evolve  from 
the  implicit  time  differencing  scheme. 

rM  _  9*h* 

=  At 

The  coarse  grid  representation  as  derived  above  are  used  to  implement  the 
multigrid  procedure.  A  flow  schematic  of  the  procedure  used  is  shown  in 
Pig.  3(b) . 

Full  weighting  is  used  for  the  f i ne-to-coarse  grid  interpolations,  taking 
into  account  all  nine  fine  grid  points  associated  with  the  coarse  grid 
equivalent  point, 


kV  -V 


Both  linear  and  third  degree  polynomial  interpolation  from  coarse  to  fine  grid 
were  tested,  and  several  different  relaxation  schemes  were  tried. 

The  Swi ten  Function 


The  values  of  the  switching  function,  g(0),  are  not  allowed  to  change 


during  a  multigrid  solution  process.  Attempts  to  let  them  vary  with  the 


solution  led  to  major  instabilities.  The  g  values  at  the  new  time  step,  gt, 


are  first  determined  from  the  previous  solution,  1 . e . ,  6  .  Using  these 

r  conv  3 


values  at  time  step  t,  the  fine  grid  equations  are  relaxed  a  certain  number  of 


times,  after  which  the  g  values  are  updated  to  the  present  9  _  .  These 

approx 


g  values  are  then  used  throughout  the  multigrid  solution. 


RESULTS 


Excellent  results  were  obtained  over  the  time  steps  prior  to  the  start  of 


cavitation,  both  in  terms  of  comparison  with  a  single  grid  iterative  solution 


and  comparison  with  Brewe's  direct  solution.  These  results  are  summarized  in 


Table  II. 


Comparisons  with  single  grid  Iteration  are  done  on  the  basis  of  work  units 


(WU)  used,  where  1  WU  is  equivalent  to  1  relaxation  sweep  over  the  finest  grid. 


Letting  M,  i.e.,  tota1  number  of  grids,  represent  the  finest  grid  and  k  a 


coarser  grid,  numbered  In  decreasing  order  respectively,  the  equivalent  WU 


used  by  gr i d  k  is 


WU^  =  4^~M) 


The  following  results  were  obtained  for  the  test  case  having  a  maximum 


eccentricity  of  0.8  and  a  minimum  eccentricity  of  0.1  (see  Table  I).  This  is 


one  of  the  more  difficult  cases  to  run,  since  very  high  pressures  and  large 
pressure  gradients  are  Induced  near  the  minimum  film  thickness.  The  pressure 
gradients  and  the  bubble  shape  change  relatively  rapidly  with  time. 


In  all  of  the  trials,  a  system  of  three  grids  was  used.  The  addition  of 
a  fourth  coarsest  grid  had  a  negligible  effect.  For  comparison  purposes,  a 
single  grid  solution  using  Gauss-Seidel  relaxation  on  a  96-  by  24-point  mesh 
was  run  (96  points  circumferentially  and  24  points  axially).  The  96-  by 
24-point  mesh  ensures  a  grid  aspect  ratio  (AR  -  Ax/Az)  of  1 .  It  required  about 
300  WU  per  time  step. 

Gauss-Seidel  (G-S)  and  Jacobi  (J)  relaxation  schemes  with  no 
overrelaxation  were  found  to  be  the  most  effective  smoothers  for  this  problem. 
Circumferential  line  relaxation,  i.e.,  solving  simultaneously  each  line  of 
points  in  the  circumferential  direction,  is  an  effective  smoother,  but  is  not 
worth  the  substantially  greater  computational  time  needed  to  solve  for  the 
periodic  boundary  conditions,  which  introduce  corner  terms  to  the  tridiagonal 
matrix.  Line  relaxation  is  used  as  a  local  smoother,  however,  when  cavitation 
develops.  Both  the  G-S  and  J  relax  points  in  the  direction  of  the  flow, 
i.e.,  the  circumferential  direction,  sweeping  across  the  axial  direction. 
Sweeping  across  the  circumferential  direction  is  not  very  effective,  nor  is  a 
combination  of  axial  and  circumferential  sweeps.  A  red-black  scheme  is  also 
not  very  effective. 

The  difference  between  the  G-S  and  the  J  schemes  when  used  in  the 
multigrid  process  is  extremely  small.  J  relaxation  uses  an  average  of  0.5  WU 
more  than  G-S  per  time  step.  The  reason  seems  to  be  that  the  J  multigrid 
uses  the  same  number  of  fine  and  medium  grid  sweeps  as  does  the  G-S  per 
solution  and  makes  up  for  its  lower  efficiency  by  using  a  greater  number  of 
the  coarsest  grid  sweeps,  which  are  very  cheap.  The  advantage  of  using  J 
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when  working  with  a  parallel  processing  computer  is  that  the  Inner  loop  Is 
vectorizable,  since  all  points  In  a  line  are  substituted  at  the  end  of  the  line 
instead  of  point  by  point.  Using  the  CRAY  XMP  on?  J  sweep  takes  one-tenth 
the  CPU  time  of  a  G-S  sweep. 

A  multigrid  solution  using  G-S  relaxation  along  the  direction  of  flow 
and  a  linear  interpolation  scheme  from  the  coarse  to  fine  grid  was  first 
tested.  The  solution  on  a  48-  by  48-point  fine  mesh  (AR  s  3.1)  required  an 
average  of  24  WU  per  time  step.  The  solution  on  the  96-  by  24-point  mesh 
required  an  average  of  14  WU  per  time  step,  which  is  nearly  22  times  faster 
than  the  single  grid  solution. 

A  third  degree  polynomial  Interpolation  scheme  from  coarse  to  fine  grid 
was  also  tested.  Using  the  same  96-  by  24-point  mesh  and  G-S,  this  scheme 
reduces  the  work  per  time  step  to  an  average  of  7.5  WU,  half  the  work  used  by 
the  linear  scheme.  Also,  the  third  degree  polynomial  routine  takes  virtually 
the  same  amount  of  CPU  time  as  the  linear  routine  on  the  CRAY  XMP,  making  it 
highly  worthwhile.  This  scheme  used  approximately  1/40  of  the  work  used  by 
single  grid  iteration. 

The  adaptive  multigrid  cycle  was  used  to  obtain  the  above  results.  To 
determine  the  efficiency,  however,  a  V-cycle  was  also  run.  Each  cycle  reduces 
the  error  by  nearly  an  order  of  magnitude. 

The  results  also  compare  well  with  Brewe's  direct  numerical  solution,  both 
in  accuracy  and  CPU  time.  Also,  load  capacities  were  compared  at  various  time 
steps.  The  greatest  difference  found  between  the  load  capacity  values  is  two 
parts  in  10  000.  Both  the  direct  and  the  multigrid  codes  are  vectorized  to  the 
highest  efficiency.  Both  were  run  on  the  CRAY  XMP  for  5000  time  steps  of 
uncavitated  flow.  The  direct  solution  took  1086  sec  CPU,  while  the  multigrid 
code  took  57  sec  CPU,  about  1/20  the  CPU  time  of  the  direct  solution. 


The  presence  of  cavltated  points  in  the  flow,  i.e.,  the  presence  of  an 
area  having  g  -  0  bounded  by  points  having  g  =  1,  requires  a  more  involved 
approach.  On  a  single  grid,  the  algorithm  handles  the  cavitation  as 
efficiently  as  it  would  an  uncavitated  configuration.  Problems  begin  to  occur 
when  coarser  grids  are  added. 

Initially,  the  coarse  grid  cavitation  area  was  determined  by  injecting 
corresponding  fine  grid  g's  directly  to  the  coarse  grid  points.  Figure  4 
shows  graphs  of  the  residual  function  at  similar  states  in  the  solution  for 
both  an  uncavitated  and  a  cavltated  configuration  using  this  scheme.  These 
graphs  were  obtained  with  no  extra  smoothing  around  the  cavltated  area.  As  can 
be  seen,  high-frequency  local-error  terms  exist  around  the  cavitated  boundary, 
whereas  the  uncavitated  region  has  already  been  well  smoothed.  If  the  program 
Is  allowed  to  continue  from  this  point,  it  moves  to  the  coarser  grid,  where  the 
cavitated  boundary  residuals  are  not  well  represented.  Depending  on  how 
unsmooth  the  boundary  residuals  are,  V-cycle  results  range  from  40-percent 
error  reduction  per  cycle  to  a  slight  divergence  of  error  terms  per  cycle. 

Extra  local  smoothing  helps  Immensely,  as  would  be  expected.  The  best 
results  are  obtained  by  using  a  local  circumferential  line  relaxation  scheme 
over  the  cavitated  region  and  boundary  points.  This  scheme  is  a  very  powerful 
smoother  and  is  also  expedient,  since,  as  a  local  smoother,  it  reduces  to  a 
purely  tridiagonal  matrix  of  a  relatively  small  number  of  points. 

The  problem  still  remains  of  deciding  how  many  local  smoothing  sweeps  is 
"enough,"  or  whether  any  are  necessary  at  all.  If  the  number  of  sweeps  is  set 
such  that  the  most  difficult  cavitation  configurations  converge  efficiently, 
then  configurations  having  smoother  initial  residuals  become  much  less 
efficient.  Some  type  of  smoothness  Indicator  is  necessary.  The  present 


routine  takes  the  fine  grid  residual  function,  fk,  and  Interpolates  It  down  to 
grid  k-1  using  the  same  full  weighting  as  in  an  actual  grid  switch.  The 
coarse  grid  values  are  then  Interpolated  linearly  back  up  to  the  fine  grid  k. 
These  then  represent  "smooth"  fine  grid  residuals  and  are  compared  with  the 
actual  residuals.  If  any  of  the  actual  terms  fall  outside  of  an  envelope 
placed  around  the  smooth  term,  the  problem  is  deemed  unsmooth  and  local 
smoothing  i s  done . 

The  above  procedure  does  much  to  stabilize  the  solution  process,  resulting 
in  an  average  usage  of  40  WU  per  time  step  solution.  For  the  g-injection 
model,  however,  order  of  magnitude  error  reduction  per  cycle  is  not  usually 
obtained,  and  certain  cavitation  configurations  do  occur  which  are  very 
difficult  or  impossible  to  solve. 

As  mentioned  earlier,  the  coarse-grid  operator  Lk_1  need  not  be  the  same 
as  the  fine-grid  operator  LM.  Because  of  this  fact,  some  latitude  in  handling 
the  values  of  the  g  coefficient  in  the  coarse  grid  equations  is  permissible. 
This  led  to  a  coarse-grid  determination  scheme  for  the  g  values  that  not  only 
circumvented  occurrences  of  instability  due  to  injection  but  had  a  major 
beneficial  effect  on  the  solution.  Recall,  the  g  values  on  the  finest  grid 
are  determined  by  the  value  of  6  at  each  nodal  point,  i.e.,  g  has  a  value  of 
1  at  full  film  points  and  a  value  of  0  at  cavitated  points.  It  was  found  that 
stability  of  the  solution  across  all  possible  cavitation  configurations  can  be 
obtained  by  defining  a  parameter  FG  as: 

FG  *  [Vi,j-i  +  gi-i,j  +  gi-i,j+i  +  gi,j-i  +  gi,j  +  gi,j+i  +  gui,j-i 


gi+l , j  +  gi+l , j  +  1. 


If 

FG  =  0; 

9k-' 

If 

FG  *  0; 

In  other  words,  a  fine  grid  point  must  have  a  value  of  0  and  must  be 
surrounded  (all  eight  points)  by  points  having  gk  =  o  in  order  for  the 
corresponding  coarse  grid  point  to  be  set  to  g^-1  =  0.  Other  schemes  were 
found  to  work  but  were  not  as  efficient.  This  scheme  resulted  in  an  average  of 
20  WU  per  time  step,  the  number  of  WU  ranging  from  11  to  35  WU.  The  solution 
process  remains  stable  throughout  bubble  formation  and  bubble  collapse.  The 
time  steps  which  required  the  most  work  units  occurred  at  the  very  beginning  of 
bubble  formation,  when  there  were  very  few  cavitated  points,  and  during  bubble 
collapse.  While  the  bubble  is  collapsing,  it  is  also  experiencing  its  greatest 
amount  of  movement  downstream,  so  it  might  be  this  movement  rather  than  the 
process  of  collapse  which  requires  more  solution  time. 

A  V-cycle  analysis  shows  that  better  than  an  order  of  magnitude  error 
reduction  per  cycle  is  obtained,  though  at  the  cost  of  extra  smoothing  on  the 
finer  grids. 

The  results  for  cavitated  flow  also  compare  well  with  Brewe's  direct 
solution.  The  same  bubble  shape,  motion,  and  duration  are  obtained  from  both 
programs.  Figure  5  contains  computed  pressure  distributions  at  various  time 
steps.  The  cavitated  area  is  indicated  by  the  flat  area  of  zero  pressure 
gradient.  Comparison  of  load  capacity  terms  shows  a  maximum  difference  of  five 
parts  in  10  000  in  the  cavitated  region.  When  run  on  the  CRAY  XMP,  the  direct 
solution  of  5000  time  steps  of  cavitated  flow  again  takes  1086  sec  CPU.  The 
multigrid  so'ution  of  5000  time  steps  of  cavitated  flow  takes  150  sec  CPU, 
approximately  one-eighth  the  CPU  taken  by  the  direct  solution.  Even  though  the 
cavitated  configurations  take  more  CPU  than  do  the  uncavitated  configurations, 
the  multigrid  solution  still  represents  a  very  significant  savings  over  the 
direct  method. 
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In  his  paper  on  the  direct  ADI  solution  of  the  Elrod  algorithm,  Brewe  (2) 
states  that  his  solution  uses  two  to  three  times  the  computational  time  used 
by  an  iterative  solution  of  a  nonconservative  film  model  using  pseudo-Gumbel 
boundary  conditions.  This  nonconservative  film  model  is  only  suited  to  steady- 
state  conditions,  but  is  often  used  by  industry.  Thus,  the  multigrid  solution 
of  the  Elrod  algorithm  requires  about  one-tenth  to  one-third  the  computational 
time  (for  uncavitated  and  cavitaced  flow  respectively)  of  the  nonconservative 
film  model  solution,  while  still  retaining  the  more  realistic  representation  of 
the  flow. 

CONCLUSION 

A  multigrid  iterative  technique  is  used  in  the  solution  of  the  Elrod 
algorithm  for  the  case  of  a  dynamically  loaded  journal  bearing  undergoing 
cavitation.  This  solution  is  compared  both  to  a  single  grid  iterative  solution 
in  terms  of  work  used,  and  to  a  direct  ADI  solution  in  terms  of  computer  time 
required.  Excellent  results  are  obtained  both  prior  to  and  during  cavitation, 
although  the  presence  of  cavitation  does  introduce  difficulties  in  the  solution 
process . 

The  best  results  are  obtained  using  the  following:  a  grid  aspect  ratio  of 
1;  full  weighting  interpolation  from  the  fine  grid  to  the  coarse  grid;  third 
degree  polynomial  Interpolation  from  the  coarse  grid  to  the  fine  grid;  either 
Gauss-Seidel  or  Jacobi  relaxation  with  no  overrelaxation.  Implementing  these 
techniques,  the  solution  at  time  steps  prior  to  cavitation  uses  1/40  the  amount 
of  work  used  by  a  single  grid  iterative  Gauss-Seidel  solution  and  1/20  the 
computer  time  used  by  the  direct  ADI  solution.  During  cavitation,  the 
multigrid  solution  uses  1/15  the  amount  of  work  used  by  a  single  grid  G-S 
solution  and  one-eighth  the  computation  time  used  by  the  direct  ADI  solution. 


Based  on  the  results  stated  In  this  paper,  it  is  evident  that  the  solution 
of  the  Elrod  algorithm  using  multigrid  techniques  provides  an  extremely  viable 
method  to  industry  for  the  solution  of  journal  bearing  problems. 

REFERENCES 

1.  Dowson,  D.  and  Taylor,  C.M.,  1979,  "Cavitation  in  Bearings,"  Annual  Review 
of  Fluid  Mechanics,  Vol .  11,  M.  VanDyke,  J.V.  Wehausen,  and  J.L.  Lumley, 
eds.,  Annual  Reviews  Inc.,  Palo  Alto,  CA,  pp.  35-66. 

2.  Brewe,  D.E.,  1986,  "Theoretical  Modeling  of  the  Vapor  Cavitation  in 
Dynamically  Loaded  Journal  Bearings,"  Journal  of  Tribology,  Vol.  108, 

No.  4,  pp.  628-638. 

3.  Jakobsson,  B.,  and  Floberg,  L.,  1957,  "The  Finite  Journal  Bearing 
Considering  Vaporization,"  Rept.  No.  190,  Chalmers  University  of 
Technology,  Goteborg,  Sweden. 

4.  Olsson,  K.O.,  1965,  "Cavitation  in  Dynamically  Loaded  Journal  Bearings," 
Rept.  No.  308,  Chalmers  University  of  Technology,  Goteborg,  Sweden. 

5.  Elrod,  H.G.  and  Adams,  M.L.,  1975,  "A  Computer  Program  for  Cavitation  and 
Starvation  Problems,"  Cavitation  and  Related  Phenomena  in  Lubrication, 

D.  Dowson,  M.  Godet,  and  C.M.  Taylor,  eds.,  Mechanical  Engineering 
Publications,  New  York,  pp.  37-42. 

6.  Elrod,  H.G.,  1981,  "A  Cavitation  Algorithm,"  Journal  of  Lubrication 
Technology,  Vol.  103,  No.  3,  pp.  350-354. 

7.  Jacobson,  B.O.,  and  Hamrock.  B.J.,  1983,  "Vapor  Cavitation  in  Dynamically 
Loaded  Journal  Bearings,"  2nd  International  Conference  on  Cavitation, 


Mechanical  Engineering  Publications,  New  York,  pp.  133-140. 

8.  Jacobson,  B.O.,  and  Hamrock,  B.J.,  1983,  "Hign-Speed  Motion  Picture  Camera 
Experiments  of  Cavitation  in  Dynamically  Loaded  Journal  Bearings,"  Journal 
of  Lubrication  Technology,  Vol.  105,  No.  3,  pp.  446-452. 


9.  Brandt,  A.,  1977,  "Mul ti -Level  Adaptive  Solutions  to  Boundary-Value 
Problems,"  Mathematics  of  Computation,  Vol .  31,  No.  138,  pp.  333-390. 

10.  Brandt,  A.,  1984,  "Multigrid  Techniques:  1984  Guide;  With  Applications 
to  Fluid  Dynamics,"  Computational  Fluid  Dynamics,  Vol.  1,  VKI-LS-1 984-04, 
Von  Karman  Institute  for  FLuid  Dynamics,  Rhode  St.  Genese,  Belgium,  1984, 
Paper  No.  1 . 

11.  lubrecht,  A. A.,  Ten  Napel  ,  W.E.,  and  Bosma,  R.,  1986,  "Multigrid,  An 
Alternative  Method  for  Calculating  Film  Thickness  and  Pressure  Profiles  in 
El astohydordynami cal ly  Lubricated  Line  Contacts,"  Journal  of  Tribology, 

Vol .  108,  No.  4,  pp.  551-556. 


TABLE  I.  -  OPERATING  CONDITIONS 


AR,  m 
R,  m 
L/D  . 


us,  rad/s 
aid,  rad/s 
3,  N/m2  . 
li,  N-s/m2 
pa,  N/m2 
pc,  N/m2 


.  5.0xl0~4 
.  .  0.0425 
.  .  .  1.0 
0.1  to  0.8 
.  .  -19.5 
.  .  -92.7 
.  1 . 72x1 09 
.  .  0.066 
1 ,0133x10s 
.  .  .  0.0 
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